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Abstract. To develop an understanding of recent experiments on the turbulence- 
induced melting of a periodic array of vortices in a thin fluid film, we perform a direct 
numerical simulation of the two-dimensional Navier-Stokes equations forced such that, 
at low Reynolds numbers, the steady state of the film is a square lattice of vortices. 
We find that, as we increase the Reynolds number, this lattice undergoes a series of 
nonequilibrium phase transitions, first to a crystal with a different reciprocal lattice and 
then to a sequence of crystals that oscillate in time. Initially the temporal oscillations 
are periodic; this periodic behaviour becomes more and more complicated, with 
increasing Reynolds number, until the film enters a spatially disordered nonequilibrium 
statistical steady that is turbulent. We study this sequence of transitions by using fluid- 
dynamics measures, such as the Okubo- Weiss parameter that distinguishes between 
vortical and extensional regions in the flow, ideas from nonlinear dynamics, e.g., 
Poincare maps, and theoretical methods that have been developed to study the melting 
of an equilibrium crystal or the freezing of a liquid and which lead to a natural set 
of order parameters for the crystalline phases and spatial autocorrelation functions 
that characterise short- and long-range order in the turbulent and crystalline phases, 
respectively. 
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1. Introduction 

A crystal melts into a disordered liquid when the temperature is raised beyond the 
melting point Tm] and when the liquid is cooled below Tm it freezes again. This melting 
or freezing transition, one of the most common equilibrium phase transitions, has been 
studied extensively; and the Ramakrishnan-Yussouff density-functional theory [U E|, EH EJ 
[5], which uses a variational free energy, has led to a good understanding of such freezing. 
It is natural to ask whether there are nonequilibrium analogues of this transition. 
One example is shear-induced melting of colloidal crystals [6]. We investigate another 
example of a nonequilibrium transition in which dynamically generated turbulence plays 
the role of temperature and disorders a crystal consisting of an array of vortices, of 
alternating sign, imposed on a thin fluid film by an external force. 

Recent experiments (TJ IE] have explored this transition from a low- Reynolds-number 
nonequilibrium vortex crystal, imposed on a thin fluid film by a force that is periodic in 
space, to a high-Reynolds-number, disordered, nonequilibrium liquid-type phase. This 
problem has been studied earlier by using linear-stability analysis [91 [10] and direct 
numerical simulations (DNS) (TTJ H2] . The former is well-suited to the study of the first 
instability of the vortex crystal with increasing Reynolds number Re; the latter have 
(a) studied the route to chaos, by using techniques from nonlinear dynamics jTU [12] to 
analyse the temporal evolution of this system, or (b) have mimicked [13] the experiments 
of Ref. [7] to study the creation and annihilation of vortices. 

We revisit the problem of turbulence-induced melting of a vortex crystal by 
conducting a direct numerical simulation (DNS). Our study combines methods from 
turbulence, nonlinear dynamics, and statistical physics to elucidate the nature of 
nonequilibrium phases and transitions in this system. We use the two-dimensional (2D) 
Navier-Stokes (NS) equation with air-drag-induced Ekman friction to model the thin 
fluid films used in experiments [?J El EH] ; this is a good model for flows in such films so 
long as the Mach number is small and the corrections arising from the finite thickness 
of the film and from the Marangoni effect can be neglected [Tol [ToT [T7] . We force this 
2D NS equation in a manner that mimics the forcing used in experiments and which 
yields, at low Re, a stationary, periodic array of vortices of alternating signs; we refer 
to this array as the vortex crystal. 

We then investigate the stability of this crystal as we increase the amplitude of 
the force and, therefore, the Reynolds number. Measures from nonlinear dynamics, 
including time-series, power spectra, and Poincare-type maps, can be used fruitfully 
here to examine the temporal behaviour of this system as it undergoes a sequence of 
transitions. This part of our study complements the work of Refs. [HI [12]. 

Furthermore, we elucidate the natures of the transitions from the spatially ordered 
crystal to the disordered, turbulent state by borrowing ideas from the density-functional 
theory of the freezing of a liquid into a crystal [H EJ [H |5] . Since the density field of a 
crystal is periodic, this theory uses the coefficients in the Fourier decomposition of the 
density as order parameters. Specifically, in a conventional crystal, p admits the Fourier 
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decomposition 

^( r ) = ^pGexp(tG • r), (1) 

G 

where the sum is over the vectors G of the reciprocal lattice; in the density-functional 
theory of freezing [U |3] the Fourier coefficients pc are identified as the order parameters 
of the liquid-to-crystal transition since their mean values vanish in the liquid phase 
for all nonzero reciprocal lattice vectors G. Moreover, pc=o is the mean density and, 
because most liquids are nearly incompressible, it undergoes a very small change at this 
transition. Spatial correlations can be characterized conveniently by the autocorrelation 
function g(r) = ([p(x + r)p(x)]), where the angular brackets denote Gibbsian thermal 
averages and the overline denotes spatial averaging over x; the Fourier transform of 
g(r) is related [U E] to the static structure factor 5(k). In an isotropic liquid phase the 
arguments of g and S are, respectively, r =| r | and k —\ k |. 

We use the Okubo- Weiss field A = det(A) as the analogue of the density p(r) in 
the density-functional theory of freezing; here A is the velocity-derivative matrix that 
has components = diUj, with Uj the j th component of the velocity For an inviscid, 
incompressible 2D fluid the local flow topology can be characterized via the Okubo- Weiss 
criterion [T8l IT9] ; this provides a useful measure of flow properties even if viscosity and 
Ekman friction are present [TU [TTJ [20]: Regions with A > and A < correspond, 
respectively, to centres and saddles [T8l [19] . Thus, in the nonequilibrium vortex crystal, 
A(r) is a periodic function; so, like p(r) in a conventional crystal, it admits the Fourier 
decomposition 

A(r) = ^A k exp(zk-r), (2) 

k 

where the sum is over the reciprocal-lattice vectors k; it is natural to think of Ak as the 
order parameters that characterise the vortex crystal. In terms of these order parameters 
we can define the analogue of the static structure factor S(k) for a conventional crystal; 
for the vortex crystal this is the two-dimensional spectrum 

E A (k) = (A k A_ k ); (3) 

here the angular brackets do not imply a Gibbsian thermal average, as in equilibrium 
melting, but denote an average over the nonequilibrium state of our system; as 
we show below, this state can be steady in time, or it can oscillate periodically 
or quasiperiodically, or it can be a chaotic state that is statistically steady. The 
autocorrelation function 

G(t) = (A(x + r)A(x)), (4) 

is related to E\(k) by a spatial Fourier transform. The turbulent phase is nearly 
isotropic so, to a good approximation, G depends only on r =| r | in this phase; and 
it characterises the short-range order in the system exactly as g(r) does in an isotropic 
liquid. 
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In addition to identifying a natural set of order parameters for the turbulence- 
induced melting of a two-dimensional vortex crystal in a thin fluid film, our study 
yields several interesting results that we describe qualitatively below: In particular, 
we find that, as we increase the Reynolds number, a rich sequence of transitions leads 
from the steady ordered vortex crystal to the disordered, turbulent state. The third 
column in Table [T] shows the types of nonequilibrium phases we encounter: There 
is SX, the original, steady square crystal imposed by the force; this is followed by 
SXA, steady crystals that are distorted, via large-scale spatial undulations, relative 
to SX; these give way to distorted crystals that oscillate in time, either periodically 
(OPXA) or quasiperiodically (OQPXA); finally the system becomes disordered and 
displays spatiotemporal chaos and turbulence (SCT). OPXA and OQPXA are actually 
a collection of several (perhaps an infinity) of nonequilibrium crystalline phases that 
oscillate in time; given the resolution of our calculation we can identify only some of 
these, as we describe in detail below. The combination of techniques that we use helps 
us to elucidate the natures of these phases in far greater detail than has been attempted 
hitherto. We believe that, by using the ideas we develop here, it should be possible to 
extend recent experiments [7] on turbulence-induced melting to uncover the rich series 
of nonequilibrium phase transitions that we have summarised in Table 1. 

The remaining part of this paper is organised as follows. In Sec. 2 we present 
the equations we use to model turbulence-induced melting in thin fluid films and the 
numerical methods we use. Section 3 contains a detailed description of our results. 
Section 4 is devoted to a discussion and to conclusions. 

2. Model and Numerical Methods 

We begin with the 2D Navier-Stokes (NS) equations that can be written, as follows, in 
the nondimensional form of Ref . [21] : 



Here u = (—d y ip, d x ip), ip, and to = V x u are, respectively, the velocity, stream function, 
and vorticity at the position x and time t; we choose the uniform density p = 1; a is 
the non-dimensionalised Ekman- friction coefficient, v is the kinematic viscosity, and 
= — n 3 [cos(nx) + cos(ra/)]/f2, is the non-dimensionalised force with injection wave 
vector n, Q = nF amp /(u 2 k 3 ), and a = nua'k/ F ampi where F amp is the forcing amplitude, 
a' is the Ekman friction, and lengths are non-dimensionalised via a factor k/n with k 
a wave vector or inverse length [21] . We denote the x and y components of the velocity 
as «i = « and u<i = v, respectively. The spatially periodic force F^ yields, at low Q, 
a vortex crystal that is also referred to as a cellular flow. A linear-stability analysis of 
this flow indicates that it has a primary instability [TU] at a critical Reynolds number 
Re c = a/2 which translates into a threshold value f2 Sjn = nRe c . This primary instability 
yields another vortex crystal, which is steady in time but whose unit cell is larger than 



(<9 f + u • VV = V 2 tu/fi + F w - au; 



(5) 
(6) 
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that of the original vortex crystal [TQl E2] . 

We solve Eqs. ^ and ^ numerically by using a pseudo-spectral method with a 2/3 
dealiasing cut-off and a second-order Runge-Kutta scheme for time marching (23] El] 
with a time step St = 0.01. We use N 2 collocation points; in most of our studies we use 
N = 128; we have checked in representative cases that our results are unchanged if we 
use N = 256. Our main goal has been to obtain long time series for several variables (see 
below) to make sure that the temporal evolution of our system is obtained accurately; 
most of our runs are at least as long as 3 x 10 6 <5t. We monitor the time-evolution of 
(a) the total kinetic energy E(t) = u 2 , (b) the stream function ip, (c) the vorticity u, 
(d) the Okubo- Weiss parameter A, and (e) the k = (1,0) component of the Fourier 
transform v of the y component v of u. Given these time series we obtain E\(k) at 
representative times and G(r), which is obtained by averaging over 20 configurations of 
A(r) separated from each other by 10 5 St, after transients in the first 10 6 time steps have 
been removed. From the time series of E(t) we obtain its temporal Fourier transform 
E(f) and thence the spectrum | E(f) | that helps us to distinguish between periodic, 
quasiperiodic, and chaotic temporal behaviours. We also augment this charaterisation 
by using Poincare-type sections in which we plot Q^no) versus ffivno) at successive times 
(see, e.g., Ref. [21] for the Kolmogorov flow). 

We show below that the vortex crystal melts, as we increase Q, via a rich sequence 
of transitions. The principal effect of the Ekman friction is to delay the onsets of these 
transitions; we have checked this explicitly in some cases. However, to make contact 
with earlier linear-stability and DNS studies of this problem, the results we present 
below have been obtained with no Ekman friction. Our qualitative conclusions are not 
affected by this. 

So long as Q < O s>n , the steady-state solution [2TJ . indicated by the subscript s, of 
Eq. ^ is LU Sjn = — n[cos(nx) + cos(m/)]. We examine the destabilisation of this state, 
with increasing Q, for two representative values of n, namely, n = 4 and n = 10 for which 
J1 S) 4 ~ 5.657 and f2 S) io — 14.142; in our DNS we choose the initial vorticity field to have 
the form to = ^ S; „ + 10~ 4 Y^ 1=0tm2= o[sin(m 1 x+m2y)+cos(m 1 x+m 2 y)}ml/ a/ [m\ + m§); 
and then we let the system evolve under the application of the force F u . We increase Q 
from f2 s n to 6.85S7 Sjri . in steps of 0.15fl, s>n , for n = 4 (runs R4-1 to R4-7 in Table[T]), and 
from Q Stn to 15.9f2 Si „ in steps of 0.1fi S) „, for n = 10 (runs R10-1 to R10-4 in Table [TJ. 
To trace the transition to chaos accurately from SXA to SCT we have also conducted 
runs where f2 is increased in steps of ~ 0.08f2 Si „ for n = 4 (runs R4-4 to R4-6 in Table [I]) 
and ~ 0.07f2 s „ for n — 10 (runs R10-2 and R10-3 in Table [TJ. We have benchmarked 
our numerical scheme by comparing results from our code with those of Ref. [21] , which 
deals with a Kolmogorov flow imposed by an external force of the form F^ = ncos(ny). 

For Q < ft s>n , the A field shows alternating centres and saddles, arranged in a 
two-dimensional square lattice, which we illustrate via pseudocolour plots for n = 4 
and n = 10 in Figs. [T] (a) and (b), respectively. [These patterns are reminiscent of a 
two-dimensional version of a perfectly ordered binary alloy, with two kinds of atoms, 
whose analogues here are centres and saddles.] We show corresponding pseudocolour 
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Table 1. Table indicating the number of the Run (e.g., R4-1), the values of n and 
f2, and the type of order. Here SX denotes the original, steady square crystal imposed 
by the force (the precise pattern depends on n) ; SXA denotes steady crystals that are 
distorted, via a large-scale undulation, relative to SX; OPXA indicates a crystal that is 
distorted slightly with respect to SX and which oscillates periodically in time; OQPXA 
is like OPXA but with quasiperiodic oscillations; SCT denotes the disordered phase 
with spatiotemporal chaos and turbulence. 

plots of ip in Figs. [2] (a) and (b), respectively. 
3. Results 

In this section we present the results of our numerical simulations for n = 4 and n — 10 
and the ranges of parameters given in Table [1} In the first subsection we present 
our results for n = 4; the next subsection contains our results for n = 10; in the 
last subsection we describe our results for the order parameters and spatial correlation 
functions. 

3.1. The case n = 4 

When we increase ft beyond ft Sjn =4, the steady-state solution u Stn= 4 becomes unstable. 
From the range of values of ft in our runs R4 — 2 (Table [T]) we observe that a new 
steady state is attained, which we illustrate, for ft = 6.5, via pseudocolour plots of if) 
and A in Figs. [3]^a) and (b), respectively. The new steady state is also a vortex crystal; 
however, it is different from the original vortex crystal as can be seen especially clearly 
by comparing the pseudocolour plots of ip in Figs. [2] (a) and|3^a). This difference also 
shows up as a very slight distortion of the crystalline structure in the pseudocolour plot 
of A shown in Fig. ^b). To use the language of solid state physics, this is an example 
of a very weak structural phase transition. Normally such a phase transition is mirrored 
in new superlattice peaks that appear in the reciprocal-space spectrum E\ in addition 
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Figure 1. (Colour online) Pseudocolour plots, illustrating the vortex crystal for 
< Cl Sjn , of the Okubo-Weiss field A for (a)n = 4, and (b) n — 10. Given our 
colour bar, vortical regions, i.e., centres, appear red whereas strain-dominated regions, 
i.e., saddles, appear dark blue. 
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Figure 2. (Colour online) Pseudocolour plots, illustrating the vortex crystal for 
< Q s ,n> of the streamfunction field i/j for (a)n = 4, and (b) n = 10. 
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to the dominant peaks associated with the original crystal structure; however, given the 
weakness of the distortion, such superlattice peaks are not visible, with our resolution 
in Fig. [3|c). Clear examples of such superlattice peaks appear as we increase Vt (see 
below) . 

At = 8.202, a new regime appears (run R4 — 3). The time-series of the energy 
E(t) now shows a periodic array of spikes. This regime has no analogue in a conventional 
crystal; it is a crystal that oscillates periodically in time and, to that extent, it can be 
thought of as a spatiotemporal crystal. The time between successive spikes is very large 
(__ 10 4 <5t) as shown by the plot of E(t) in Fig. |4[a); this is why our DNS runs must be 
very long to distinguish such states from one that is steady; we have also checked that 
the time between successive spikes is the same (to three-figure accuracy) for iV = 64 and 
iV = 128. In Figs. [4](b) and (c) we show pseudocolour plots of ip and A, respectively; the 
former shows a large-scale spatial undulation and the latter some deformation relative 
to the original crystal. This deformation is also mirrored in the distortion, relative 
to Fig. |3](c), of the dominant peaks in the reciprocal-space spectrum E\(k) shown in 
Fig. id). 

For runs R4 — 4, i.e., 9.05 < Q < 15.3, we find a new crystalline state that is steady 
in time. It has a large-scale spatial undulation relative to the original vortex crystal 
as illustrated, for = 11.3, by the pseudocolour plots of ip and A in Figs. [5^a) and 
(b), respectively. This undulation leads to a distortion of the dominant peaks in the 
reciprocal-space spectrum E\(k) of Fig. [5^c), which also shows new superlattice peaks 
that occur at smaller values of (k x ,k y ) relative to the dominant peaks. 

On further increasing Q we enter a new regime (runs R4 — 5, i.e., 15.3 < Q < 17.3) 
in which we have a spatiotemporal crystal, i.e., a spatially periodic A that oscillates in 
time. The time-series E(t) displays a periodic array of spikes as shown in Fig. |6^a); this 
leads to the frequency-space spectrum | E(f) | of Fig. [6^b). The peaks in this spectrum 
can be labelled as Efo, where £ is a positive integer and fo is the fundamental frequency 
that can be obtained from the inverse of the temporal separation between successive 
spikes in Fig. [6^a); this is a clear signature of periodic temporal evolution. The Poincare- 
type map in the (3f^[-0(l, 0)], 3 [f)(l, 0)]) plane, Fig. |6](c), shows that the projection of the 
attractor on this plane is a closed loop in this case. Pseudocolour plots of ip and A 
[ Figs. [7](a)-(b)] are similar, respectively, to those in Figs. [5^a)-(b) if we look at their 
spatial patterns; however, they oscillate in time as can be seen most clearly from their 
animated versions [mpeg files psi_movie_R5.mpeg-lam_movie_R5.mpeg]. The associated 
reciprocal-space spectrum __a also oscillates between the spectra shown in Figs. (7^c) 
and (d) as can be seen clearly from its animated version [avi file lamf_movie_R5.avi]. 

The time between successive spikes in E (t) decreases as Q increases. To quantify 
this, we define the inter-spike interval Tj as follows: Tj = £ i+1 — 1{, where t$ is the time 
at which E(t) crosses, for the i th time, its mean value, (E(t)), from below; we can think 
of i as the spike index. In Fig. Ma) we plot Tj versus i; this shows that the mean value 
of Tj decreases as _2 increases [Fig. [8^b)]; furthermore, Tj oscillates slightly about its 
mean value for any given value of Q. The magnitude of these oscillations, which we have 
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Figure 3. (Colour online) Pseudocolour plots for — 6.5 of (a) the streamfunction 
ip and (b) the Okubo- Weiss parameter A with superimposed contour lines, (c) A filled 
contour plot of the reciprocal-space energy spectrum E\ showing clear, dominant peaks 
at the forcing wave vectors. 
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Figure 4. (Colour online) Plot for = 8.202 of (a) the time evolution of the energy 
E(t) for — 8.202. Pseudocolour plots of (b) the streamfunction ip and (c) the Okubo- 
Weiss parameter A with superimposed contour lines, (d) A filled contour plot of the 
reciprocal-space energy spectrum E\ showing the distortions of peaks at the forcing 
wave vectors. 

used for the error bars in Fig. |8^b), decreases as fl increases. We have checked explicitly 
that our results here do not change when we increase the resolution of our DNS from 
iV = 128 to N = 256. 

At Q — 17.8, i.e., run R4 — 6, another transition occurs: The time series of the 
energy E(t) and its frequency spectrum | E(f) \ are shown, respectively, in Figs.^a) and 
(b). The latter displays peaks superimposed on a noisy background signal; these peaks 
can be indexed as f , /i, f\ — 2f , f — 2/i, and 3f — 2/i, within our numerical accuracy 
and with fo ~ 0.001653 and f\ ~ 0.001707. Since fo/fi is not a simple rational number, 
we conclude that Fig. |9|b) indicates principally quasiperiodic temporal evolution with 
a small chaotic admixture, the former associated with the peaks indexed above and the 



Turbulence-induced melting of a nonequilibrium vortex crystal 



12 




Figure 5. (Colour online) Pseudocolour plots for = 11.3 of (a) the streamfunction 
tp and (b) the Okubo-Weiss parameter A with superimposed contour lines, (c) A 
filled contour plot of the reciprocal-space energy spectrum E\ showing clear, dominant 
peaks, at the forcing wave vectors, and subdominant supcrlattice peaks at smaller wave 
vectors. 
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Figure 6. (Colour online) Plots for = 15.3 of (a) the time evolution of the energy 
E(t), (b) | E(f) | versus the frequency /, and (c) the Poincare-type section in the 
($Ru[l,0],St>[l,0]) plane. 
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latter with the noisy background signal. We believe the chaotic part of the signal comes 
from transitions between the elliptical islands in the Poincare-type section of Fig. [9|c). 
The plot of the inter-spike interval T« versus the spike index i in Fig. [9](d) confirms the 
complicated temporal evolution of this state. 

As we increase Q further (R4 — 7) the temporal evolution of the system becomes 
ever more chaotic; this is associated with a disordered pattern of vortices in space 
too. Thus we obtain a state with spatiotemporal chaos and turbulence, which is our 
analogue of the liquid state. We illustrate this for Q = 50. We begin with the time 
series of E(t) and the spectrum | E(f) | in Figs. [lO^a) and (b), respectively; the latter 
shows clearly a broad background that is indicative of temporal chaos. This is further 
confirmed by the nearly uniform spread of points in the Poincare-type section [Fig. [To|c)] 
in the (3?[w(l, 0)], ^[^(l, 0)]) plane. The disordered spatial organisation of this state is 
illustrated by the pseudocolour plots of if) and A [Figs. 11 a)-(b), respectively] and the 
reciprocal-space spectrum -E'A(k) of Fig. [TT^c) that shows several new modes in addition 
to the original peaks, whose vestiges are still visible since F w continues to act on the 
system. 

Thus we see that the turbulence-induced melting of our nonequilibrium vortex 
crystal is far richer than its equilibrium counterpart. For the case n = 4 investigated 
above it proceeds as described in column 4 in Table [TJ Before we present an analysis 
of the disordered state in terms of the spatial autocorrelation function G(r) and the 
evolution of the order parameters (Ak) with Q, we give below a short summary of our 
results for n = 10; the route to turbulence is slightly different for this case than for 
n = 4. 



3.2. The case n=10 

Our results for n = 10 are based on the runs R10 — 1 to RIO — 5 in Table [TJ 

For Q < Q s n the steady vortex crystal is shown by the pseudocolour plot of A 
in Fig. [ijb). As we increase f2 beyond f2 Si „ we find in runs R10 — 2, i.e., the range 
Q St n < Q < 22.6, a new steady state in which pseudocolour plots of if) and A show large- 
scale spatial undulations caused by small deformations of the original vortex crystal 
[Figs. [l2] (a) and (b), respectively]; consequently the dominant peaks in the reciprocal- 
space spectrum E\ are slightly distorted. 

Around Q = 24 (run RIO — 3) another transition occurs: the time series of E(t) 
is periodic [Fig. |l3^a)] and its spectrum | E(f) | [Fig. [l3|(b)] shows one dominant 
peak, i.e., higher harmonics are nearly absent. Thus the Poincare-type plot in the 
(3fcy[l, 0], 0]) plane [Fig. [l3jc)] displays a simple attractor. The spatial structure 
of this state is illustrated by the pseudocolour plots of if> and A shown, respectively, in 
Figs. [l4|a) and (b); the associated reciprocal-space spectrum Ea is shown in Fig. [l4^c). 
Given the temporal behaviour of this state, these structures, in real or reciprocal space, 
oscillate in time at the frequency given by the temporal evolution of E(t). For Q = 24 
the large-scale spatial structures in if) oscillate around their mean positions; but, for 
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25 < Vt < 28, we find a travelling- wave type pattern, which reenters our simulation 
domain by virtue of the periodic boundary conditions that we use in our pseudo-spectral 
method. If we compare the frequency spectra | E{f) | for the cases Q = 24 and Q = 28, 
we find higher harmonics in the latter but they are all multiples of one fundamental 
frequency. 



In Fig. 15 a) we show plots of Tj versus i (cf. Fig. uB[a) for n = 4) for various values 



of VL. From these we obtain the plot of the mean value (Tj) versus VL shown in Fig. 15 
This first decreases, as we increase Q, and then increases mildly at Q = 28. 

For Q > 29 the time-series of E(t) appears chaotic and the associated frequency 



spectrum | E(f) | displays a broad background, as we show in the illustrative Figs. 16 



(a) and (b) for f2 = 225. The associated Poincare-type section in Fig. 16 c) confirms 
that the temporal behaviour is chaotic. The spatial patterns are also disordered as we 
show via the pseudocolour plots of ip an d A in Figs. [I7|(a) and (b), respectively. The 
corresponding reciprocal-space spectrum E\ shows that a large number of modes are 
excited. Thus we have both spatial disorder and temporal chaos; as in the case n = 4, 
the analogue of the liquid state is a turbulent one with spatiotemporal chaos. However, 
given the resolution in Q that we have been able to obtain in our calculations, the route 
to this state of spatiotemporal chaos is different for n = 10 and n = 4 as can be seen 
from column 4 in Table Q] 

3. 3. Order parameters and spatial autocorrelation functions 

We now return to ideas borrowed from the density-functional theory [TJ [3J, HJ [5] of 
freezing by examining the behaviour of the order parameters (A^) as functions of Q. 
Recall that, in equilibrium melting, pc jumps discontinuously from a nonzero value in 
the crystal to zero in the liquid at the first-order melting transition. As we have noted 
above, the turbulence-induced melting of our vortex crystal is far more complicated; it 
proceeds via a sequence of transitions. In turbulence-induced melting of a vortex crystal 
5f(Ak), obtained by summing A^ over the four forcing wave vectors, is the equivalent of 
5J(Pg) m a conventional crystal; 3?(Ak) changes with Q as shown, respectively, for (a) 
n = 4 and k = (4, 4) and (b) n = 10 and k = (10, 10) in Figs. [18] (a) and (b). For small 
values of Q the steady state is SX so, in Fourier space, only modes with the forcing wave 
vector are significant. On increasing Q, the crystal undergoes a series of transitions that 
take it from the crystal SX to the disordered, turbulent state SCT; as this happens the 
spectral weight slowly shifts away the forcing wave numbers. This explains the trend 
observed in Figs. [18J 

The spatial correlations in the crystalline and disorderd states can be characterised 
by the spatial autocorrelation function G(r) defined in Eq. (4). Representative plots are 



shown, respectively, for n = 4 and n = 10 in Figs. 19 (a)-(d). For the crystalline case 



we evaluate G(r) along the line connecting r = (7r/2,7r/2) and r = (7r/2, 7r); this shows 



a periodic array of peaks [Figs. 19 (a) and (b) for n = 4 and n = 10, respectively]; the 



widths of these peaks are related to the widths of vortical or strain-dominated regions. 
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In the turbulent phase we present data obtained by a circular average of G [Figs. 19 
(c) and (d) for n = 4 and n = 10, respectively]; here the peaks decay over a length 
scale that indicates the degree of short-range order. This decay is similar to the decay 
of spatial correlation functions in a disordered liquid. 

4. Conclusions 

We have carried out a detailed numerical study of turbulence-induced melting of a 
nonequilibrium vortex crystal in a forced, thin fluid film. We use ideas from the 
density-functional theory of freezing [TJ |3j HJ E], nonlinear dynamics, and turbulence 
to characterise this. Correlation functions, similar to those in liquid-state theory, have 
been used by some recent experiments |8] to analyse the short-range order in the 
turbulent phase; nonlinear-dynamics methods, such as Poincare-type maps, have been 
used in the numerical studies of Ref . [H] ; experimental studies have used the curvature 
of Lagrangian trajectories to identify extrema in vortical and strain-dominated regimes. 
To the best of our knowledge, there is no study that brings together a variety of methods, 
as we do, to analyse turbulence-induced melting. 

The advantages of our approach are as follows: (a) it helps us to identify the 
order parameters for turbulence-induced melting and thus contrast it with conventional 
melting; (b) the sequence of transitions can be characterised completely in terms of the 
Eulerian fields ip and A, the total energy E(t), and suitable Fourier transforms of these; 
(c) spatial correlations in crystalline and turbulent phases can be studied conveniently 
in terms of G. 

Equilibrium phase transitions occur strictly only in the thermodynamic limit that 
is, roughly speaking, the limit of infinite size [25J. It is interesting to ask how we might 
take the thermodynamic limit for the vortex crystals we have studied here. There seem 
to be at least two ways to do this: (a) in the first the system size can be taken to infinity 
in such a way that the areal density of the vortical and strain-dominated regimes remains 
the same in the ordered, crystalline phase; (b) in the second way we can increase the 
parameter n in the forcing F u so that more and more unit cells occupy the simulation 
domain (see, e.g., Figs.[l|a) and (b) for n = 4 and n = 10, respectively). Such issues have 
not been addressed in detail by any study, partly because, for large system sizes, it is 
not possible to obtain the long time series that are required to characterise the temporal 
evolution of the system (especially in the states we have referred to as spatiotemporal 
crystals). In particular, it is quite challenging to investigate the system-size dependence 
of the transitions summarised for n = 4 and n = 10 in Table [TJ From Figs. 19 (c) and 



(d) we can extract a correlation length; this length is much smaller than the linear size 
of our simulation domain so we expect that our results in the turbulent phase (SCT) 
will not change if we increase the size of our system. However, subtle size dependence 
might occur in the ordered phase as follows: as we increase Q, the original crystal is 
distorted by large-scale spatial undulations that are associated with the inverse cascade 
of energy in two-dimensional turbulence; if these undulations lead to crystalline phases 
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with larger and larger unit cells, then our finite-size calculations will become unreliable 
when the size of the unit cell becomes comparable to the size of our simulation domain. A 
systematic study of such subtle finite-size effects is a very challenging task that requires 
much more detailed simulations than have been attempted so far. 

As we have shown above, the sequence of transitions that comprise turbulence- 
induced melting of a vortex crystal is far richer than conventional equilibrium melting. 
There is another important way in which the former differs from the latter: To maintain 
the steady states, statistical or strictly steady, we always have a force F^, thus, in 
the language of phase transitions, we always have a symmetry-breaking field, both in 
the ordered and disordered phases. Strictly speaking, therefore, there is no symmetry 
difference between the disordered, turbulent state and the ordered vortex crystal, as 
can be seen directly from the remnants of the dominant peaks in the reciprocal-space 



spectra E\(k) in Figs. 11 and[l7[c) for n = 4 and n = 10, respectively. One consequence 
of this is that the order parameters (Ak), with k equal to the forcing wave vectors, do 
not vanish identically in the disordered, turbulent phase; however, they do assume very 
small values. Moreover, in the case of turbulence-induced melting the crystal undergoes 
a transition from an ordered state to an undulating crystal to a fully turbulent state. 
Thus there is no noise and hence no fluctuations in the crystalline state; i.e., it is 
equivalent to a crystal at zero temperature. This has no analogue in the equilibrium 
melting of a crystal. 

In equilibrium, different ensembles are equivalent; we can, e.g., use either the 
canonical or the grand-canonical ensemble to study the statistical mechanics of a system 
and, in particular, the phase transitions in it. However, this equivalence cannot be 
taken for granted when we consider nonequilibrium statistical steady states (see, e.g., 
Ref. [26]). We have seen an example of this in Ref. [TT] where certain PDFs show 
slightly different behaviours depending on whether we keep the Grashof number (i.e., 
the nondimensionalised amplitude of the force) fixed or whether we keep the Reynolds 
number fixed. Turbulence-induced melting offers another example of the inequivalence 
of dynamical ensembles: the precise sequence of transitions that we encounter in going 
from the vortex crystal to the turbulent state depends on whether we do so by changing 
the Grashof number as in Ref. [TT] or whether we do so by changing Q as we have done 
here. We have checked explicitly that we can reproduce the sequence of transitions in 
Ref. [TT] if we tune the Grashof number rather than Q to obtain the turbulent state. 

Investigations of similar transitions, such as in the Kolmogorov flow [211 127] . can 
benefit by using the combination of methods we have used above. Detailed studies of the 
effects of confinement, air-drag induced Ekman friction on turbulence-induced melting, 
initiated, e.g., in Refs. [12], [22], can also make use of our methods, but that lies beyond 
the scope of this paper. We hope, therefore, that our study will encourage experimental 
groups to analyse turbulence-induced melting by using the set of techniques and ideas 
that we have described above. 
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Figure 7. (Colour online) Pseudocolour plots for = 15.3 of (a) the streamfunction 
i\) and (b) the Okubo- Weiss parameter A with superimposed contour lines, (c) and (d) 
Filled contour plots of showing the two spectra between which it oscillates in time. 



Turbulence-induced melting of a nonequilibrium vortex crystal 



20 



3500 
3000 
2500 
2000 
1500 



(a) 



vAAAAA/WWWWWlAArtAAAAAAAAA/' 



20 40 60 80 100 120140 160 




E-T 2500 



2000 



1500 



15.5 16.0 16.5 

Q 



17.0 



Figure 8. (Colour online) Plots of (a) the inter-spike interval Tj versus the spike index 
i for il — 15.3 (red curve), = 15.8 (black curve), = 16.3 (purple curve), = 16.8 
(green curve), and £1 = 17.3 (blue curve) and (b) the time mean value of Ti versus fl 
(see text for error bars). 
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Figure 9. (Colour online) Plots for £1 — 17.8 of (a) the time evolution of the energy 
E(t), (b) the spectrum E(f) versus /, (c) the Poincare-type section in the plane 
(5R[t)(l,0)],3 i [u(l,0)]), and (d) the inter-spike interval Tj versus the spike index i. 
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Figure 10. (Colour online) Plots for = 50 of (a) the time evolution of the energy 
E(t), (b) the spectrum | E(f) | versus /, and (c) the Poincare-type section in the 
(3?[u(l,0)],3[w(l,0)]) plane. 
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Figure 11. (Colour online) Pseudocolour plots for f2 = 50 of (a) the streamfunction 
ip and (b) the Okubo- Weiss parameter A with superimposed contour lines, (c) A filled 
contour plot of the reciprocal-space energy spectrum E\ which shows that a large 
number of modes are excited. 
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Figure 12. (Colour online) Pseudocolour plots for £1 = 22.62 of (a) the streamfunction 
ip and (b) the Okubo- Weiss parameter A with superimposed contour lines, (c) A filled 
contour plot of the reciprocal-space energy spectrum E\ showing clear, but slightly 
distorted, dominant peaks at the forcing wave vectors. 
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Figure 13. (Colour online) Plots for O = 24 of (a) the time evolution of energy E(t) 
versus t, (b) the power spectrum E(f) versus /, and (c) the Poincare-type section 
S[w(l,0)] versus SJ[u(l,0)]. 
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Figure 14. (Colour online) Pseudocolour plots for f2 = 24 of (a) the streamfunction 
i\) and (b) the Okubo- Weiss parameter A with superimposed contour lines, (c) A filled 
contour plot of the reciprocal-space energy spectrum E\ showing clear, dominant peaks 
at the forcing wave vectors. 
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Figure 15. (Colour online) Plots of (a) the inter-spike interval Tj versus the spike 
index i for S7 = 24 (blue curve), SI = 25 (purple curve), SI = 26 (green curve), SI = 27 
(red curve), and SI = 28 (cyan curve) and (b) the time mean value Tj for different 
values of S7. 
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Figure 16. (Colour online) Plots for = 225 of (a) the time evolution of E(t) 
, (b) the spectrum | E(f) | versus /, and (c) the Poincare-type section in the 
(3?[i)(l,0)], $[{*(!, 0)]) plane. 
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Figure 17. (Colour online) Pseudocolour plots for = 225 of (a) the streamfunction 
ip and (b) the Okubo- Weiss parameter A with superimposed contour lines, (c) A filled 
contour plot of the reciprocal-space energy spectrum showing that a large number 
of modes are excited. 




Figure 18. (Colour online) Plot showing a decrease in (Ak) with increasing f2 for (a) 
n = 4 and k = (4,4) and (b) n = 10 and k = (10, 10). 
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Figure 19. (Colour online) Plots of G(r): (i) crystalline state: (a) n = 4,f2 < il s n 
and (b) n — 10, < fi s „ along the line connecting r = (ir/2,ir/2) and r = (ir/2,ir); 
and (ii) circularly averaged in the turbulent state: (c) n = 4, fi = 20.81 and (d) 
n = 10, ft = 225. 



